##########################################################################################
# Title: Build ALL_CITIES (2000 & 2010) — Section 2 Aggregation
# Section: 2 (Inputs → Processed_BD_Analysis)
# Purpose: Collapse per-city distance matrices into banded aggregates and save
# IO:
#   IN : Processed_Data/Intermediate/Section_2/distance_matrix/<City>_<comb>.rds
#   OUT: Processed_Data/Intermediate/Section_2/Processed_BD_Analysis/ALL_CITIES_<comb>.rds
##########################################################################################

rm(list = ls()); options(stringsAsFactors = FALSE)

suppressPackageStartupMessages({
  library(data.table)
})

# ---- Project root (align with run_all.R) -------------------------------------
if (!requireNamespace("here", quietly = TRUE)) install.packages("here", quiet = TRUE)
ROOT   <- normalizePath(here::here(), winslash = "/")
setwd(ROOT)

IN_DIR  <- file.path(ROOT, "Processed_Data", "Intermediate", "Section_2", "distance_matrix")
OUT_DIR <- file.path(ROOT, "Processed_Data", "Intermediate", "Section_2", "Processed_BD_Analysis")
dir.create(OUT_DIR, recursive = TRUE, showWarnings = FALSE)

# ---- Target combinations (ONLY these) ----------------------------------------
combinations <- c("100K_09_2000", "100K_09_2010")

# ---- Cities ------------------------------------------------------------------
cities <- c(
  "Johannesburg","Capetown","Tshwane","Durban","NM_Bay",
  "Nairobi","Lagos","Abidjan","Cairo","London","Istanbul","Seoul",
  "SaoPaulo","KualaLumpur","Moscow","Mexico_City","Buenos_Aires",
  "Shenzhen","Frankfurt","Zurich","Bogota","Tehran","Dubai",
  "New_York","Mumbai"
)

# ---- Helper: process one city into aggregated distance bands -----------------
process_city_file <- function(city, comb) {
  infile <- file.path(IN_DIR, paste0(city, "_", comb, ".rds"))
  if (!file.exists(infile)) {
    message("⚠️ Missing: ", infile)
    return(NULL)
  }
  dt <- readRDS(infile) |> as.data.table()
  
  # Expect: x, y, building, pop, ntl, min_dist_cbd, bd_dummy, building_cutoff, cluster_size, year, city
  dt <- dt[!is.na(ntl) & !is.na(pop) & pop >= 10]
  
  # Distance bands
  dt[, bounds := fcase(
    min_dist_cbd <  5000,  5,
    min_dist_cbd < 10000, 10,
    min_dist_cbd < 15000, 15,
    min_dist_cbd < 20000, 20,
    min_dist_cbd < 25000, 25,
    default = 30
  )]
  
  dt[, counts := .N, by = bounds]
  
  # Collapse by bands
  agg <- dt[, .(
    pop_sum = sum(pop, na.rm = TRUE),
    ntl_sum = sum(ntl, na.rm = TRUE),
    counts  = max(counts, na.rm = TRUE)
  ), by = bounds]
  
  agg[, pop_density := pop_sum / counts]
  agg[, ntl_density := ntl_sum / counts]
  agg[, pop_km_density := NA_real_]
  
  # References within 5 km
  pop5  <- agg[bounds == 5,  pop_sum]
  popd5 <- agg[bounds == 5,  pop_density]
  ntld5 <- agg[bounds == 5,  ntl_density]
  
  agg[, population_5    := ifelse(length(pop5),  max(pop5,  na.rm = TRUE), NA_real_)]
  agg[, rel_pop_5       := pop_sum / population_5]
  agg[, population_d_5  := ifelse(length(popd5), max(popd5, na.rm = TRUE), NA_real_)]
  agg[, rel_pop_d_5     := pop_density / population_d_5]
  agg[, rel_pop_diff_5  := (pop_density - population_d_5) * 100 / population_d_5]
  agg[, ntl_d_5         := ifelse(length(ntld5), max(ntld5, na.rm = TRUE), NA_real_)]
  agg[, rel_ntl_d_5     := ntl_density / ntl_d_5]
  agg[, rel_ntl_diff_5  := (ntl_density - ntl_d_5) * 100 / ntl_d_5]
  agg[, city_name := city]
  
  agg
}

# ---- Main: build ONLY the two ALL_CITIES files -------------------------------
for (comb in combinations) {
  
  # Phase 1: per-city aggregation (in memory)
  per_city_list <- lapply(cities, process_city_file, comb = comb)
  per_city_list <- per_city_list[!vapply(per_city_list, is.null, logical(1))]
  if (!length(per_city_list)) {
    message("No cities processed for ", comb)
    next
  }
  
  # Phase 2: combine cities
  all_dt <- rbindlist(per_city_list, use.names = TRUE, fill = TRUE)
  
  # City-level population at each band (b = 5,10,15,20,25)
  for (b in seq(5, 25, 5)) {
    all_dt[bounds == b, paste0("pop_b", b) := pop_sum]
    all_dt[, paste0("population_b", b) := max(get(paste0("pop_b", b))), by = city_name]
  }
  
  # Flags
  all_dt[, south_africa := city_name %in% c("Capetown","Johannesburg","Durban","Tshwane","NM_Bay")]
  all_dt[, other_africa := city_name %in% c("Lagos","Nairobi","Abidjan","Cairo")]
  
  # Display names (paper)
  rename_map <- c(
    "Capetown"      = "Cape Town",
    "Mexico_City"   = "Mexico City",
    "SaoPaulo"      = "São Paulo",
    "KualaLumpur"   = "Kuala Lumpur",
    "Tshwane"       = "Tshwane (Pretoria)",
    "Durban"        = "eThekwini (Durban)",
    "NM_Bay"        = "NM Bay (Port Elizabeth)",
    "Buenos_Aires"  = "Buenos Aires",
    "New_York"      = "New York City"
  )
  all_dt[, city_name := ifelse(city_name %in% names(rename_map), rename_map[city_name], city_name)]
  
  # Global means outside SA & Other Africa
  glob <- all_dt[!south_africa & !other_africa,
                 .(global_rel_ntl_diff_5 = mean(rel_ntl_diff_5, na.rm = TRUE),
                   global_rel_pop_diff_5 = mean(rel_pop_diff_5, na.rm = TRUE)),
                 by = bounds]
  all_dt <- merge(all_dt, glob, by = "bounds", all.x = TRUE)
  
  # Derived fields
  all_dt[, diff_global_ntl    := rel_ntl_diff_5 - global_rel_ntl_diff_5]
  all_dt[, diff_global_pop    := rel_pop_diff_5 - global_rel_pop_diff_5]
  all_dt[, pop_km_density     := pop_density * 100]  # GHSL 100m → per km²
  all_dt[, ln_pop_km_density  := log(pmax(pop_km_density, .Machine$double.eps))]
  all_dt[, city_code          := as.integer(factor(city_name))]
  all_dt[, bounds_f           := relevel(factor(bounds), ref = "20")]
  
  # Phase 3: save ONLY the ALL_CITIES file
  out_file <- file.path(OUT_DIR, paste0("ALL_CITIES_", comb, ".rds"))
  saveRDS(all_dt, out_file)
  message("✅ Wrote: ", normalizePath(out_file, winslash = "/"))
}

message("Done: 100K_09_2000 and 100K_09_2010 ALL_CITIES files.")

